# Space separator
# brew install pandoc
knitr::opts_chunk$set(warning = FALSE) # setting of rmarkdown
options(width = 100)
# render('heart_integration.Rmd')
sc_raw <- read.table("../asp_heart/count/Filtered/scRNA/all_cells_count_matrix_filtered.tsv",
sep = "\t", as.is = T)
sc_meta <- read.table("../asp_heart/count/Filtered/scRNA/all_cells_meta_data_filtered.tsv",
sep = "\t", as.is = T, fill = T, header = T) # pending for meta-data problem
sc_mx <- data.matrix(sc_raw[-1, -1]) # this is raw matrix
colnames(sc_mx) <- sc_raw[1, -1]
rownames(sc_mx) <- sc_raw[-1, 1]
sc_tp <- sapply(colnames(sc_mx), function(x) {
sc_meta[sc_meta[, 1] == x, "celltype"]
})
sum(is.na(sc_tp)) #[1] 0
# Xu's data
load("../../2022/result/filter2022/clustering/hvgBysam500_mg2_noNao1_tf291_tfLcc_hvgLcc_rmEndog3/flc_raw_0421.rdata") # THIS IS RAW matrix of ALL USED CELLS
load(file = "../../2022/result/ficlu2022/ficlu_20220522.rdata")
in_tp <- c("splanchnic LPM-12", paste("heart", 1:11, sep = "-"), "endothelium-4",
"endothelium-6", "endothelium-10") # include endocardium
length(vis_blood <- rownames(allc_meta)[allc_meta[, "dissection_part"] %in%
c("viscera", "viscera-lowerTrunk") & allc_meta[, 3] == "blood"]) # 3023
our_mx <- flc_raw[, c(rownames(allc_anno)[allc_anno[, 1] %in% in_tp], vis_blood)] # [1] 32351 8266
rm(flc_raw)
gene <- read.table("../../10x/head/filtered_gene_bc_matrices/filtered_gene_bc_matrices/GRCh38/genes.tsv",
sep = "\t", as.is = T) # [1] 33694 2
ge2an <- gene[, 2]
names(ge2an) <- gene[, 1]
get_id <- function(x) {
sapply(x, function(y) {
names(ge2an)[ge2an == y][1]
})
}
# merge matrix
get_emb <- function(x) {
emb <- gsub("[A-z]", "", sapply(x, function(y) {
strsplit(y, split = "_")[[1]][1]
}))
return(emb)
}
# remove blood cells in Asp's data
length(asp_rm <- names(sc_tp)[sc_tp %in% c("Erythrocytes", "Immune cells")]) # DO NOT remove
# Heather's data
load("../etchevers/code/et_data.rdata")
sum(!comg %in% rownames(et_raw)) # [1] 0, checked on genes!
# normalize matrix for 3 datasets
our_norm <- t(t(our_mx) * 10000/colSums(our_mx))
asp_norm <- t(t(sc_mx) * 10000/colSums(sc_mx))
etc_norm <- t(t(et_raw) * 10000/colSums(et_raw))
save(our_norm, asp_norm, etc_norm, file = "h3_norm_mxs.rdata")
length(comg <- names(ge2an)[names(ge2an) %in% rownames(our_mx) & ge2an %in%
rownames(sc_mx)]) # [1] 15108
mg_mx <- list(our_mx[comg, !sapply(colnames(our_mx), get_emb) %in% c("21",
"22")], our_mx[comg, sapply(colnames(our_mx), get_emb) %in% c("21",
"22")], sc_mx[ge2an[comg], ], et_raw[comg, et_meta[, 1]]) # separate our embryos into two batches
rownames(mg_mx[[3]]) <- rownames(mg_mx[[1]])
sapply(mg_mx, dim)
# [2,] 6829 1437 3777 49227
rm(et_raw)
rm(sc_raw, sc_mx)
rm(our_mx)
# meta data
allc_meta <- read.table(file = "../../GSE157329/allc_meta_order.txt", sep = "\t",
header = T, as.is = T)
rownames(allc_meta) <- allc_meta[, 1]
our_meta <- cbind("our", allc_meta[colnames(our_mx), c(2, 5, 8, 3)])
colnames(our_meta)[c(1, 5)] <- c("dataset", "system")
our_meta[our_meta[, 5] != "blood", 5] <- "heart"
asp_meta <- cbind("asp", sc_tp[colnames(sc_mx)], sc_tp[colnames(sc_mx)],
"6.5w", "heart")
colnames(our_meta) -> colnames(asp_meta)
asp_meta[asp_meta[, 3] %in% c("Erythrocytes", "Immune cells"), 5] <- "blood"
# add heather's meta data
et_st <- c("8.6w", "9.0w", "10.7w")
names(et_st) <- c("hu140_8.6_pcw", "hu084_9.0_pcw", "hu122_10.7_pcw")
et_meta2 <- cbind("etc", as.character(et_meta[, c("celltype")]), as.character(et_meta[,
c("celltype")]), as.character(et_st[et_meta[, "stage"]]), "heart")
et_meta2[et_meta2[, 2] %in% c("Blood"), 5] <- "blood"
rownames(et_meta2) <- et_meta[, 1]
colnames(et_meta2) <- colnames(our_meta)
# all meta data
all_meta <- rbind(our_meta, asp_meta, et_meta2)
sum(!unlist(sapply(mg_mx, colnames)) %in% rownames(all_meta)) # 0, check
######################## UPDATE CNC and immune cell labels in Asp
######################## data, which is swapped!
all_meta[all_meta[, 1] == "asp" & all_meta[, 3] == "Immune cells", 3] <- "token"
all_meta[all_meta[, 1] == "asp" & all_meta[, 3] == "Cardiac neural crest cells ",
3] <- "Immune cells"
all_meta[all_meta[, 1] == "asp" & all_meta[, 3] == "token", 3] <- "Cardiac neural crest cells "
vis_rmg <- scan("../../cross_embryo/code/vis_rmg.txt", what = character(0))
tf <- scan("../../Imai_table/humanTF_GO0003700_clear.txt", what = character(0))
# Seurat
library(Seurat)
batch <- "blood_hvg1k" # the best
seu <- lapply(mg_mx, function(x) {
CreateSeuratObject(x, min.cells = 0, min.features = 0)
})
seu <- mapply(function(x, y) {
if (y)
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 1000)
}, x = seu, y = c(T, T, T, T))
length(features <- SelectIntegrationFeatures(object.list = seu, nfeatures = 1000))
length(features <- setdiff(features, vis_rmg))
seu <- lapply(X = seu, FUN = function(x) {
x <- ScaleData(x, features = features, verbose = FALSE)
x <- RunPCA(x, features = features, verbose = FALSE)
})
anchor <- FindIntegrationAnchors(object.list = seu, anchor.features = features,
reduction = "rpca", dims = 1:50, reference = 4)
seu_comb <- IntegrateData(anchorset = anchor)
DefaultAssay(seu_comb) <- "integrated"
# Run the standard workflow for visualization and clustering
seu_comb <- ScaleData(seu_comb, verbose = FALSE)
seu_comb <- RunPCA(seu_comb, npcs = 50, verbose = FALSE, features = features)
seu_comb <- RunUMAP(seu_comb, reduction = "pca", dims = 1:30)
umap <- seu_comb@reductions$umap@cell.embeddings
pcac <- seu_comb@reductions$pca@cell.embeddings
save(pcac, all_meta, tp2_cncc, file = paste("../etchevers/result/", batch,
"_pcac.rdata", sep = ""))
# clustering
seu_comb <- FindNeighbors(seu_comb, dims = 1:30)
seu_clu <- rep(0, nrow(umap))
resos <- seq(0.1, 1, 0.1)
for (reso in resos) {
seu_comb <- FindClusters(seu_comb, resolution = reso, dims.use = 1:30)
one_res <- as.numeric(seu_comb@active.ident)
names(one_res) <- names(seu_comb@active.ident)
seu_clu <- cbind(seu_clu, one_res)
}
seu_clu <- seu_clu[, -1]
colnames(seu_clu) <- resos
# calculate markers for level-1 clusters at resolution = .1
reso <- 0.1
seu_comb <- FindClusters(seu_comb, resolution = reso, dims.use = 1:30)
seu_mk <- FindAllMarkers(seu_comb, only.pos = T)
rm(seu_comb, seu, anchor)
# plot seurat clusters
# pdf(paste('../etchevers/result/',batch,'_seuclu.pdf',sep=''),
# width=8,height=12)
par(cex = 2, las = 1, mar = c(5, 4, 1, 1), lwd = 3, mfrow = c(2, 2), pch = 16)
for (i in 1:8) {
plot_cen <- sapply(sort(unique(seu_clu[, i])), function(x) {
cl <- rownames(seu_clu)[seu_clu[, i] == x]
res <- c(median(umap[cl, 1]), median(umap[cl, 2]))
return(res)
})
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1",
ylab = "umap-2", frame = F, col = col6[-1][seu_clu[, i]][rand_ind],
cex = 0.2, main = max(seu_clu[, i]))
for (i in 1:ncol(plot_cen)) text(plot_cen[1, i], plot_cen[2, i], label = i,
xpd = NA)
}
# dev.off()
.
# plot by stage and cell types
tmp <- unique(read.table("/Volumes/baolab-1/bao_data_zrc/baolab/xuy/ciona/cluster2/product/code/type_all_col3.txt",
sep = "\t", as.is = T, comment.char = "")[, 2])
col6 <- c("lightgrey", "red", tmp[3], "purple", tmp[2], tmp[5], "blue",
"green4", "orchid", "turquoise", "sienna", tmp[9], "yellow2", "hotpink",
"navy", "steelblue", "skyblue", "pink", "black", tmp[4], rainbow(7))
st_col <- col6[c(18, 20, 2, 8, 17, 7, 9)]
names(st_col) <- c("CS12", "CS13-14", "CS15-16", "6.5w", "8.6w", "9.0w",
"10.7w")
sys_col <- col6[c(3, 2)]
names(sys_col) <- c("heart", "blood")
tp1_col <- col6[c(2, 4, 9, 14, 27, 18, 20, 7, 16, 17, 6, 5, 3, 23, 8)]
names(tp1_col) <- unique(all_meta[all_meta[, 1] == "our" & all_meta[, 5] !=
"blood", 3])
tp1_col <- tp1_col[c(1:7, 11, 8:10, 13:15, 12)]
# tp1_col[15]<-NA
tp2_col <- col6[c(1, 16, 10, 5, 8, 11, 5, 3, 7, 14, 9, 19)]
names(tp2_col) <- unique(all_meta[all_meta[, 1] == "asp", 3])
tp2_col <- tp2_col[c(11, 10, 6, 9, 2, 3, 8, 5, 4, 12)]
tp3_col <- c(tp1_col[c(2, 2, 15, 4, 4, 4, 8, 3, 3, 8, 12, 14, 9, 10, 10,
10, 10, 13, 8)], "black", "grey")
names(tp3_col) <- sort(unique(all_meta[all_meta[, 1] == "etc", 3]))
rand_ind <- sample.int(nrow(umap))
# pdf(paste('../etchevers/result/',batch,'_stage.pdf',sep=''),width=7,height=7)
# # by dataset
par(cex = 2, las = 1, mar = rep(0.5, 4), lwd = 3)
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1", ylab = "umap-2",
frame = F, col = st_col[as.character(all_meta[rownames(umap), 4])][rand_ind],
cex = 0.2, xaxt = "n", yaxt = "n")
legend("bottomright", col = st_col, legend = c(names(st_col)[1:3], "6.5-7w",
names(st_col)[5:7]), xpd = NA, bty = "n", pch = 16, cex = 0.8)
# dev.off()
.
# pdf(paste('../etchevers/result/',batch,'_system.pdf',sep=''),width=7,height=7)
# # by system
par(cex = 2, las = 1, mar = rep(0.5, 4), lwd = 3)
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1", ylab = "umap-2",
frame = F, col = sys_col[as.character(all_meta[rownames(umap), 5])][rand_ind],
cex = 0.2, xaxt = "n", yaxt = "n")
legend("bottomright", col = sys_col, legend = names(sys_col), xpd = NA,
bty = "n", pch = 16, cex = 1.5)
# dev.off()
.
plot_cen <- sapply(names(tp1_col), function(x) {
cl <- rownames(all_meta)[all_meta[, 3] == x]
res <- c(median(umap[cl, 1]), median(umap[cl, 2]))
return(res)
})
plot_cen[2, 13] <- plot_cen[2, 13] - 1.2
plot_cen[1, 8] <- plot_cen[1, 8] - 0.5
plot_cen[2, 1] <- plot_cen[2, 1] - 0.5
# pdf(paste('../etchevers/result/',batch,'_tp1.pdf',sep=''),width=7,height=7)
par(cex = 2, las = 1, mar = rep(0.5, 4), lwd = 3)
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1", ylab = "umap-2",
frame = F, col = tp1_col[as.character(all_meta[rownames(umap), 3])][rand_ind],
cex = 0.2, xaxt = "n", yaxt = "n")
for (i in 1:15) text(plot_cen[1, i], plot_cen[2, i], label = i, xpd = NA)
legend("bottomright", col = tp1_col, legend = paste(1:15, c("second heart field",
names(tp1_col)[-c(1, 15)], "splanchnic LPM-21")), xpd = NA, bty = "n",
pch = 16, cex = 0.4)
# dev.off()
.
plot_cen <- sapply(names(tp2_col), function(x) {
cl <- rownames(all_meta)[all_meta[, 3] == x]
res <- c(median(umap[cl, 1]), median(umap[cl, 2]))
return(res)
})
# pdf(paste('../etchevers/result/',batch,'_tp2.pdf',sep=''))
par(cex = 2, las = 1, mar = rep(0.5, 4), lwd = 3)
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1", ylab = "umap-2",
frame = F, col = tp2_col[as.character(all_meta[rownames(umap), 3])][rand_ind],
cex = 0.2, xaxt = "n", yaxt = "n", main = "")
for (i in 1:ncol(plot_cen)) text(plot_cen[1, i], plot_cen[2, i], label = i,
xpd = NA)
legend("bottomleft", col = tp2_col, legend = paste(1:length(tp2_col), names(tp2_col)),
xpd = NA, bty = "n", pch = 16, cex = 0.4)
# dev.off()
.
sum(duplicated(unlist(sapply(unique(all_meta[, 1]), function(x) {
unique(all_meta[all_meta[, 1] %in% x, 3])
})))) # [1] 0, no duplication on cell type name across dataset
## [1] 0
plot_cen <- sapply(names(tp3_col)[-3], function(x) {
# not include blood
cl <- rownames(all_meta)[all_meta[, 3] == x]
res <- c(median(umap[cl, 1]), median(umap[cl, 2]))
return(res)
})
# pdf(paste('../etchevers/result/',batch,'_tp3.pdf',sep=''))
par(cex = 2, las = 1, mar = rep(0.5, 4), lwd = 3)
plot(umap[rand_ind, 1], umap[rand_ind, 2], pch = 16, xlab = "umap-1", ylab = "umap-2",
frame = F, col = tp3_col[as.character(all_meta[rownames(umap), 3])][rand_ind],
cex = 0.2, xaxt = "n", yaxt = "n", main = "")
for (i in 1:ncol(plot_cen)) text(plot_cen[1, i], plot_cen[2, i], label = i,
xpd = NA)
legend("bottomright", col = tp3_col[-3], legend = paste(1:(length(tp3_col) -
1), names(tp3_col)[-3]), xpd = NA, bty = "n", pch = 16, cex = 0.4)
# dev.off()
.
# batch<- 'blood_hvg1k'
batch <- "blood_hvg1k_cncc2" # separate Asp's CNCC as two clusters
all_meta[names(tp2_cncc), 3] <- tp2_cncc
# Slingshot
library(slingshot)
tp <- unname(all_meta[rownames(pcac), 3])
table(tp)
slin <- getLineages(data.matrix(pcac[, 1:30]), tp) # ~30 mins
# write.table( slin@slingParams$dist
# ,file=paste('../etchevers/',batch,
# '_sling_dist.txt',sep=''),quote=F,sep='\t')
slin_dist <- read.table("../etchevers/result/blood_hvg1k_cncc2_sling_dist.txt",
sep = "\t", as.is = T)
rownames(slin_dist) -> colnames(slin_dist)
slin_dist <- data.matrix(slin_dist)
cl_gp[[4]][[2]] <- c("cncc-1", "cncc-2") # update 2 cncc clusters in Asp's data
# Slingshot distance and convert to similarity
plot_mx <- 1/(1 + slin_dist[unlist(cl_gp), unlist(cl_gp)]) # 1/(1+dist): best
plot_max <- 0.25
plot_mx[plot_mx > plot_max] <- plot_max
# plot_min<- 2; plot_mx[plot_mx< -plot_min]<- plot_min
diag(plot_mx) <- plot_max
plot_lab <- rownames(plot_mx)
plot_lab[plot_lab == "cardiomyocyte-2"] <- "splanchnic LPM-21"
plot_lab[plot_lab == "Secondary heart field (SHF)"] <- "second heart field"
# pdf(paste('../etchevers/result/',batch,'_all_tp_slin_dist.pdf',sep='')
# )
tmp <- heatmap3(plot_mx, labRow = plot_lab, labRow2 = NA, scale = "none",
dendrogram = "none", trace = "none", Rowv = F, Colv = F, symkey = F,
density.info = "none", keysize = 1, col = colorRampPalette(c("blue",
"white", "red"))(499), color_key_label = "similarity", color_key_label_cex = 1,
margins = c(1, 1), color_key_axis_cex = 1, key_mar = c(2, 0.5, 0.5,
0.5), labCol = NA, labRow_pos = c(2, 4), sepwidth = c(0.1, 0.1),
sepcolor = "black", cexRow = 0.7, cexCol = 1, lhei = c(1, 5), lwid = c(1,
5), labCol_las = 1, labCol_pos = 3, RowSideColors = ds_col[tp_ds[rownames(plot_mx)]],
ColSideColors = ds_col[tp_ds[rownames(plot_mx)]], colsep = c(18, 31,
40), rowsep = c(18, 31, 40))
# dev.off()
.
# 1) Best match and second best match z-score > XX; All match require
# > YY; 2) Only match Etc to Asp and Asp to our data
slin_sim <- 1/(1 + slin_dist[unlist(cl_gp), unlist(cl_gp)])
slin_zs <- lapply(rownames(slin_sim), function(x) {
vec <- slin_sim[x, setdiff(colnames(slin_sim), x)]
(vec - mean(vec))/sd(vec)
}) # exclude self dist
names(slin_zs) <- rownames(slin_sim)
zs_cut <- 1
slin_cut <- 0.06
slin_link <- sapply(names(slin_zs), function(x) {
if (tp_ds[x] == "our")
return(NULL)
if (tp_ds[x] == "asp")
tp <- names(slin_zs[[x]])[tp_ds[names(slin_zs[[x]])] == "our"] else tp <- names(slin_zs[[x]])[tp_ds[names(slin_zs[[x]])] == "asp"]
res <- NULL
best <- tp[which.max(slin_zs[[x]][tp])]
if (slin_sim[x, best] > slin_cut)
res <- best
if (slin_zs[[x]][tp[order(slin_zs[[x]][tp], decreasing = T)][2]] >
zs_cut) {
next_best <- tp[order(slin_zs[[x]][tp], decreasing = T)][2]
if (slin_sim[x, next_best] > slin_cut)
res <- c(res, next_best)
}
return(res)
})
slin_link[["Ventricular cardiomyocytes"]] <- c(slin_link[["Ventricular cardiomyocytes"]],
"ventricle cardiomyocyte-1")
slin_hit <- matrix(0, nr = nrow(slin_sim), nc = ncol(slin_sim))
rownames(slin_hit) <- rownames(slin_sim) -> colnames(slin_hit)
for (i in 1:length(slin_link)) {
if (length(slin_link[[i]]) > 0) {
slin_hit[names(slin_link)[i], slin_link[[i]]] <- 1
slin_hit[slin_link[[i]], names(slin_link)[i]] <- 1
}
}
plot_mx <- slin_hit
# pdf(paste('../etchevers/result/',batch,'_all_tp_slin_dist_cut.pdf',sep='')
# )
tmp <- heatmap3(plot_mx, labRow = rownames(plot_mx), labRow2 = NA, scale = "none",
dendrogram = "none", trace = "none", Rowv = F, Colv = F, symkey = F,
density.info = "none", keysize = 1, col = colorRampPalette(c("blue",
"white", "red"))(499), color_key_label = "hit", color_key_label_cex = 1,
margins = c(1, 1), color_key_axis_cex = 1, key_mar = c(2, 0.5, 0.5,
0.5), labCol = NA, labRow_pos = c(2, 4), sepwidth = c(0.1, 0.1),
sepcolor = "black", cexRow = 0.7, cexCol = 1, lhei = c(1, 5), lwid = c(1,
5), labCol_las = 1, labCol_pos = 3, RowSideColors = ds_col[tp_ds[rownames(plot_mx)]],
ColSideColors = ds_col[tp_ds[rownames(plot_mx)]], colsep = c(18, 31,
40), rowsep = c(18, 31, 40))
# dev.off()
.
# output hit by igraph
library(igraph)
ds_sp <- ds_col
ds_sp[1:3] <- c("circle", "triangle", "square")
# 4 groups
cl_gp_mx <- matrix(0, nr = 4, nc = length(unlist(cl_gp)))
rownames(cl_gp_mx) <- c("CM", "Epi", "Endo", "CNCC")
colnames(cl_gp_mx) <- unlist(cl_gp)
for (i in 1:length(cl_gp)) cl_gp_mx[i, unlist(cl_gp[[i]])] <- 1
hit_4ig <- cbind(rbind(matrix(0, nr = 4, nc = 4), t(cl_gp_mx)), rbind(cl_gp_mx,
slin_hit))
rownames(hit_4ig)[1:4] <- colnames(hit_4ig)[1:4]
hit <- graph_from_adjacency_matrix(hit_4ig, mode = "undirected")
V(hit)$color <- c(rep(NA, 4), ds_col[tp_ds[rownames(slin_hit)]])
V(hit)$frame.color <- c(rep(NA, 4), ds_col[tp_ds[rownames(slin_hit)]])
V(hit)$label.cex <- 0.5
V(hit)$label.color <- "black"
V(hit)$label <- c(rep("", 4), rownames(slin_hit))
V(hit)$label[V(hit)$label == "cardiomyocyte-2"] <- "splanchnic LPM-21"
V(hit)$label[V(hit)$label == "Secondary heart field (SHF)"] <- "second heart field"
V(hit)$shape <- c(rep("circle", 4), ds_sp[tp_ds[rownames(slin_hit)]])
E(hit)$color <- "black"
E(hit)[.inc(c("CM", "Epi", "Endo", "CNCC"))]$color <- NA
for (i in 1:length(slin_link)) {
if (length(slin_link[[i]]) > 0)
E(hit)[slin_link[[i]][1] %--% names(slin_link)[i]]$color <- "red"
}
E(hit)$width <- 3
hit_lay = layout_nicely(hit)
# pdf(paste('../etchevers/result/',batch,'_all_tp_slin_dist_igraph.pdf',sep=''),
# width=6,height=6 ) par(mar=c(.5,2,.5,.5)) plot(hit, vertex.size=8 ,
# vertex.label.dist=1, vertex.label.degree =pi,layout=hit_lay,
# vertex.label.cex=.7 ) legend('topleft', col=ds_col,
# legend=c('CS12-CS16','6.5-7w','8.6-10.7w'),xpd=NA,pch=c(16,17,15) )
# legend('topright', col=c('red','black'), legend=c('best match',
# '2nd match'),xpd=NA,lty=1,lwd=3) dev.off()
# triangle shape https://igraph.org/r/doc/shapes.html
mytriangle <- function(coords, v = NULL, params) {
vertex.color <- params("vertex", "color")
if (length(vertex.color) != 1 && !is.null(v)) {
vertex.color <- vertex.color[v]
}
vertex.size <- 1/150 * params("vertex", "size")
if (length(vertex.size) != 1 && !is.null(v)) {
vertex.size <- vertex.size[v]
}
symbols(x = coords[, 1], y = coords[, 2], bg = vertex.color, fg = "NA",
stars = cbind(vertex.size, vertex.size, vertex.size), add = TRUE,
inches = FALSE)
}
# clips as a circle
add_shape("triangle", clip = shapes("circle")$clip, plot = mytriangle)
# new layout as tree
cl_gp_mx <- matrix(0, nr = 4, nc = length(unlist(cl_gp)))
rownames(cl_gp_mx) <- c("CM", "Epi", "Endo", "CNCC")
colnames(cl_gp_mx) <- unlist(cl_gp)
for (i in 1:length(cl_gp)) cl_gp_mx[i, cl_gp[[i]][[1]]] <- 1
hit_tree <- cbind(rbind(matrix(0, nr = 4, nc = 4), t(cl_gp_mx)), rbind(cl_gp_mx,
slin_hit))
rownames(hit_tree)[1:4] <- colnames(hit_tree)[1:4]
fk_link <- c("Mcph2", "cncc-1", "cncc-2", "Smooth muscle cells ")
names(fk_link) <- c("Endothelium / pericytes ", "cardiomyocyte-1", "Secondary heart field (SHF)",
"epicardium")
for (i in 1:length(fk_link)) {
hit_tree[names(fk_link)[i], fk_link[i]] <- 1
hit_tree[fk_link[i], names(fk_link)[i]] <- 1
}
hit_tree <- hit_tree[c(1:4, 6:5, 7, 9:8, 10:nrow(hit_tree)), c(1:4, 6:5,
7, 9:8, 10:nrow(hit_tree))] # change node order. NOTICE the order is different with 'slin_hit'
hit2 <- graph_from_adjacency_matrix(hit_tree, mode = "undirected")
V(hit2)$color <- c(rep(NA, 4), ds_col[tp_ds[rownames(slin_hit)]])
V(hit2)$frame.color <- c(rep(NA, 4), ds_col[tp_ds[rownames(slin_hit)]])
V(hit2)$label.cex <- 0.5
V(hit2)$label.color <- "black"
V(hit2)$label <- c(rep("", 4), rownames(hit_tree)[-(1:4)])
V(hit2)$label[V(hit2)$label == "cardiomyocyte-2"] <- "splanchnic LPM-21"
V(hit2)$label[V(hit2)$label == "Secondary heart field (SHF)"] <- "second heart field"
V(hit2)$shape <- c(rep("circle", 4), ds_sp[tp_ds[rownames(slin_hit)]])
E(hit2)$color <- "grey"
E(hit2)[.inc(c("CM", "Epi", "Endo", "CNCC"))]$color <- NA
for (i in 1:length(slin_link)) {
if (length(slin_link[[i]]) > 0)
E(hit2)[slin_link[[i]][1] %--% names(slin_link)[i]]$color <- "red"
}
for (i in 1:length(fk_link)) E(hit2)[fk_link[[i]][1] %--% names(fk_link)[i]]$color <- NA
E(hit2)$width <- 2
ds_pos <- c(pi, pi, 0)
names(ds_pos) <- names(ds_col)
ds_dist <- c(2, 1, 1)
names(ds_dist) <- names(ds_col)
lab_pos <- c(rep(0, 4), ds_pos[tp_ds[rownames(hit_tree)[-(1:4)]]])
lab_dist <- c(rep(0, 4), ds_dist[tp_ds[rownames(hit_tree)[-(1:4)]]])
hit_lay2 = layout_as_tree(hit2, root = 1:4)
hit_lay2 <- -hit_lay2[, 2:1]
# pdf(paste('../etchevers/result/',batch,'_all_tp_slin_dist_igraph2.pdf',sep=''),
# width=6,height=6 )
par(mar = c(0.5, 0.5, 0.5, 3))
plot(hit2, vertex.size = 8, vertex.label.dist = lab_dist, vertex.label.degree = lab_pos,
layout = hit_lay2, vertex.label.cex = 0.8)
# legend('topleft', col=ds_col,
# legend=c('CS12-CS16','6.5-7w','8.6-10.7w'),xpd=NA,pch=c(16,17,15),
# cex=.8 )
legend("bottomleft", col = c("red", "grey"), legend = c("best match", "2nd match"),
xpd = NA, lty = 1, lwd = 3, cex = 0.8)
# dev.off()
# save.image('asp_scrna_blood.rdata.rdata') # save workspace
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.43 SeuratObject_4.1.3 Seurat_4.3.0.1
## [4] slingshot_2.8.0 TrajectoryUtils_1.8.0 SingleCellExperiment_1.22.0
## [7] SummarizedExperiment_1.30.2 Biobase_2.60.0 GenomicRanges_1.52.0
## [10] GenomeInfoDb_1.36.1 IRanges_2.34.1 S4Vectors_0.38.1
## [13] BiocGenerics_0.46.0 MatrixGenerics_1.12.2 matrixStats_1.0.0
## [16] princurve_2.1.6 igraph_1.5.0 rmarkdown_2.23
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_1.8.7 magrittr_2.0.3
## [4] spatstat.utils_3.0-3 zlibbioc_1.46.0 vctrs_0.6.3
## [7] ROCR_1.0-11 spatstat.explore_3.2-1 RCurl_1.98-1.12
## [10] htmltools_0.5.5 S4Arrays_1.0.4 sass_0.4.6
## [13] sctransform_0.3.5 parallelly_1.36.0 KernSmooth_2.23-21
## [16] bslib_0.5.0 htmlwidgets_1.6.2 ica_1.0-3
## [19] plyr_1.8.8 plotly_4.10.2 zoo_1.8-12
## [22] cachem_1.0.8 mime_0.12 lifecycle_1.0.3
## [25] pkgconfig_2.0.3 Matrix_1.5-4.1 R6_2.5.1
## [28] fastmap_1.1.1 GenomeInfoDbData_1.2.10 fitdistrplus_1.1-11
## [31] future_1.33.0 shiny_1.7.4.1 digest_0.6.33
## [34] colorspace_2.1-0 patchwork_1.1.2 tensor_1.5
## [37] irlba_2.3.5.1 progressr_0.13.0 fansi_1.0.4
## [40] spatstat.sparse_3.0-2 polyclip_1.10-4 httr_1.4.6
## [43] abind_1.4-5 compiler_4.3.1 highr_0.10
## [46] MASS_7.3-60 DelayedArray_0.26.6 tools_4.3.1
## [49] lmtest_0.9-40 httpuv_1.6.11 future.apply_1.11.0
## [52] goftest_1.2-3 glue_1.6.2 nlme_3.1-162
## [55] promises_1.2.0.1 grid_4.3.1 Rtsne_0.16
## [58] cluster_2.1.4 reshape2_1.4.4 generics_0.1.3
## [61] gtable_0.3.3 spatstat.data_3.0-1 tidyr_1.3.0
## [64] data.table_1.14.8 sp_2.0-0 utf8_1.2.3
## [67] XVector_0.40.0 spatstat.geom_3.2-2 RcppAnnoy_0.0.21
## [70] ggrepel_0.9.3 RANN_2.6.1 pillar_1.9.0
## [73] stringr_1.5.0 later_1.3.1 splines_4.3.1
## [76] dplyr_1.1.2 lattice_0.21-8 deldir_1.0-9
## [79] survival_3.5-5 tidyselect_1.2.0 miniUI_0.1.1.1
## [82] pbapply_1.7-2 gridExtra_2.3 scattermore_1.2
## [85] xfun_0.39 stringi_1.7.12 lazyeval_0.2.2
## [88] yaml_2.3.7 evaluate_0.21 codetools_0.2-19
## [91] tibble_3.2.1 BiocManager_1.30.21.2 cli_3.6.1
## [94] uwot_0.1.16 xtable_1.8-4 reticulate_1.30
## [97] munsell_0.5.0 jquerylib_0.1.4 Rcpp_1.0.11
## [100] spatstat.random_3.1-5 globals_0.16.2 png_0.1-8
## [103] parallel_4.3.1 ellipsis_0.3.2 ggplot2_3.4.2
## [106] bitops_1.0-7 listenv_0.9.0 viridisLite_0.4.2
## [109] scales_1.2.1 ggridges_0.5.4 leiden_0.4.3
## [112] purrr_1.0.1 crayon_1.5.2 rlang_1.1.1
## [115] cowplot_1.1.1 formatR_1.14